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Abstract 

Biological systems (among others) may respond to 
a large variety of distinct external stimuli, or sig- 
nals. These perturbations will generally be presented 
to the system not singly, but in various combina- 
tions, so that a proper understanding of the system 
response requires assessment of the degree to which 
the effects of one signal modulate the effects of an- 
other. This paper develops a pair of structural met- 
rics for sparse differential equation models of complex 
dynamic systems and demonstrates that said met- 
rics correlate with proxies of the susceptibility of one 
signal-response to be altered in the context of a sec- 
ond signal. One of these metrics may be interpreted 
as a normalized arc density in the neighborhood of 
certain influential nodes; this metric appears to corre- 
late with increased independence of signal response. 

1 Introduction 

Biological signaling pathways frequently intersect one 
another, leading to the phenomenon of cross-talk, 
wherein the effect of one signaling pathway influ- 
ences the activity of another. For example, during 
Drosophila development, the epidermal growth fac- 
tor receptor (EGFR) and Notch signaling pathways 
interact both within and between cells, producing 
highly context-specific responses to their respective 
signals during the formation of spatially structured 
organs including the eye pQ. The interactions be- 
tween these pathways range from antagonistic to co- 
operative across different processes and at different 
times within a given process, resulting in an intricate 
integration of signal response. 

The concept of signal interaction embodied in this 
example can be generalized to less obvious biologi- 
cal examples as well: for instance, it has been sug- 
gested that cell-to-cell variation between embryonic 
stem (ES) cells may result in differential responses to 



differentiation signals, thereby permitting some ES 
cells to undergo lineage specification while others re- 
main pluripotent [2]. Regarding the differences in 
molecular population and configuration constituting 
such ES cellular variation as random perturbations 
(or more abstractly, signals) to the biochemical state, 
this phenomenon may also be described as the mod- 
ulation by one signal of the systemic response to an- 
other. 

A key question in such cases is to what extent the 
presence and magnitude of one signal interferes with 
or reinforces the effects of the other. Perhaps the sim- 
plest possibility is linear superposition, in which the 
system essentially responds with the sum of the re- 
sponses it would have to the two signals if presented 
separately — i.e., despite making use of common sig- 
naling components, the two signals act in a funda- 
mentally independent manner. Such "signal indepen- 
dence" is, however, inconsistent with many properties 
ascribed to the interlinked networks of biological sig- 
naling pathways. For instance, while an AND gate of 
sorts could be constructed with linear superposition 
of two signals if the effect of either one alone was 
below threshold, in the presence of noisy signals of 
widely varying magnitudes (more the rule than the 
exception in biological signaling) , a particularly high 
magnitude single signal would lead to activation by 
itself. Meanwhile, even the complexity of a simple 
XOR (cxclusive-or) gate would be impossible with 
linear signal superposition/signal independence. 

Previous work by Wylie [3j suggests that the net- 
work topology of dynamic systems plays a key role in 
determining the integration of multistable dynamic 
"switches:" specifically, sparse networks of relatively 
homogenous node degree were found to be more fa- 
vorable to switch integration than were dense or 
scale-free networks. Here we investigate whether sim- 
ilar ideas might be applied to wider class of signal in- 
tegration phenomena not necessarily involving multi- 
stability. 



Pertubed and Permuted 






Figure 1: Characteristic network structures: rdg stands for random digraph, sf indicates scale-free, sw and latt indicate 
small- world digraphs with rewiring probabilities p rw = 0.1 and p rw = 0, respectively. Numerical postfixes indicate average node 
in- and out-degree; where not indicated, average degree is 6. Arc directionality suppressed for visual clarity. 



2 Terminology and Notation 

We consider deterministic nonlinear dynamic systems 
with steady state at the origin, so that 

f = /,M (i) 

= A% 3 X 3 + ^2 ^'./' + 0(X 3 ) 

For the purposes of this paper, we will restrict our- 
selves to quadratic systems for which all third-order 
and higher terms vanish, so that the linearization ma- 
trix A and the three-index array B (which, without 
loss of generality, we assume is symmetric with re- 
spect to its 2 nd and 3 rd indices) totally determine the 
dynamics. We will also write equation (JlJ as 

f =Ax + B(x,x) (2) 

Models of large biological systems are generally 
sparse in the sense that most entries Aij in the matrix 
A vanish. We can thus associate a network structure 



(more precisely, a directed graph in which "one-loop" 
arcs from a node i to itself are allowed) G — (Vg, Eg) 
with the system by 

(i j) G E G iff A jt ^ (Note index order) (3) 

We will assume here that the quadratic array B is 
consistent with the network structure G, in the sense 
that 

B ljk = unless (j ->• i), {k i) e E G (4) 

which is necessary if the linearized system structure 
G is to be stable to small perturbations of the form 

f (x) h> f (x) + Ac (5) 

as is discussed further in section |4j The network 
structures G considered in this work are generally 
random digraphs, scale-free digraphs characterized 
by high variance of node (in- and out-) degree, 
and small- world digraphs (including lattice digraphs) 
characterized by high clustering. (Here in-/out- 
degree arc defined ignoring both one-loop arcs and 
arc weights.) Figure [l] offers visualizations of some 
characteristic structures. 
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Figure 2: Arc set 
K=(123)(45)(6), with k 
lighted in dark red. 



representation of 
= 6 and c(K) = 



permutation 
3, is high- 



3 Topological Properties of 
Characteristic Polynomial 

We here refer to the characteristic polynomial asso- 
ciated with the linearization matrix A, 



det(A- A7) = (-1) 



(6) 



fc=0 



as the characteristic polynomial of the system de- 
scribed by equation ([I]). The notation Fk used in 
equation (|6| is intended to evoke the interpretation 
of the coefficients of the characteristic polynomial as 
the "feedback at length fc" in the (weighted) network 
G. This interpretation stems from the relationship 

El El 



F k = 



E 

Kee k 



E wk 

Kee k 



(7) 

where is the set of all permutations of k inte- 
gers chosen from the set {1,2, ... ,n}. Here we regard 
any particular permutation K chosen from Qk as a 
set of arcs (i —¥ j) corresponding to the mappings 
of individual elements by the permutation operation 
(illustrated in figure [2|. Considering the standard 
cycle representation of permutation groups [5], it is 



apparent that, considered graphically, K will gener- 
ally consist of a definite number c(K) of node-disjoint 
cycles, the sum of whose lengths is k. We define here 
the weight wr of the permutation K with respect to 
the system A to be the product of the arc weights Aji 
associated with the arcs (i —> j) making up K (times 
the sign factor (— l) c ^' +1 ), so that Fk is simply the 
sum of the length-fc permutation cycle weights. This 
is the basis for considering Fk as a measure of feed- 
back (of length k) in the linearized system A. 

4 Variation under System Per- 
turbation 

We take signal inputs to our dynamic systems in the 
form of perturbations to the system dynamics 



f (x) H. f (x) + Ac 



(8) 



In the context of a biochemical model, signals de- 
scribed by equation (|8| might consist of the steady in- 
put and/or removal of a given set of chemical species. 

For small perturbations, the root Ay (= AyW + 
Ay( 2 -* + higher order terms) of the function (f + Ac) 
is given to first order in Ac by 



Ay 



(i) 



-A _1 Ac 



(9) 



If the system is perturbed by two distinct pertur- 
bations Ac and 5c — representing here two incom- 
ing signals — the first order root shift response 
(AyW + 8y^) will obviously be a linear superpo- 
sition of the responses to the two individual signals 
input singly. We then consider also the second order 
mixed terms, 



ASy 



(2) 



-2A- x B(Ay« *y {1) ) 



(10) 



As the only terms we will be interested in of greater 
than first order will be mixed terms of the AS 
form, we will henceforth drop the parenthetical su- 
perscripts. 

We here consider two distinct outputs, in the form 
of changes to the steady-state properties of the dy- 
namic system under consideration, resulting from the 
signal inputs Ac and Sc. These are the coefficients 
of the characteristic polynomial Fk (for brevity, we 
will generally refer simply to "the characteristic poly- 
nomial Fk") and the eigenvalues A (particularly the 
"least stable" eigenvalue (s) Ai s of largest real part 
Re(A)). These particular outputs were chosen for 
their close relationship with system stability, as well 
as for some degree of analytic convenience. 
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One simple metric for quantifying the degree to 
which a signal Ac applied to the system f may be ex- 
pected to change the effects of a randomly distributed 
second signal 8c with regard to a particular output 
function U is the Pearson correlation 



Corr 5 (SU + ASU, SU) 

1 , ((SUASU)) 2 S «[A«/] 2 )), 



(11) 



2(([6U]*)) S 2 (d SU } 2 )) S 



H.O.T. 



where the correlation and (co)variances are all taken 
with respect to the distribution of Sc (the notation 
((XY)) indicates the covariance of the random vari- 
ables X and Y). Equation (11) is a simple lowest 
order expansion of the Pearson correlation in ASU. 
It is worth noting that since we are considering a 
quantity which clearly depends on the overall scale of 
the perturbations SU, AU, and ASU, and since we 
would like to compare this quantity across systems 
with very different structures, we must consider how 
to properly normalize the signals Sc and Ac. For 
the purposes of this work, signals were normalized 
by the root-mean-square average magnitudes of the 
first-order shifts <5Ai s to the least-stable (largest real- 
part) eigenvalues of the linearization matrix A (see 
sections [6] - [7]) . 



To apply equation (11) to our cases of interest 
(U = F k or U = Ais), we must thus obtain estimates 
of SF k , ASFk, SX\ S , and A(5Ai s . These estimates are 
primarily obtained through numerical methods, but 
it is of interest to derive formulae for a few of them in 
terms of the system parameters A and B. Linearizing 
around the shifted steady state 5y, we find that the 
linearization matrix A is shifted to A + SA, with 5 A 
given by 



and 



SA, 



ASA, 



2^> 



ijk 



2^B ijk A8y k 



(12) 



(13) 



Equation ( 12 ) also provides the rationale for the con- 



sistency constraints on the entries of array B men- 
tioned in section [2] above, since for any arc (j — > 
i) £ E G , SAij = for all perturbations Sy requires 
Bijk = B ik] = for all fc. 

The shift in the system steady state to Sy under 
the perturbation Sc then induces a shift in the char- 
acteristic polynomial given by 



(14) 



while under the combination of perturbations Sy and 
Ay, 



ASF k = 



f>F k 



ASA, 



d 2 F k 
dAadA„ 



AAj i jSA qr 



(15) 

where the range of the second summation in equation 



(15) excludes the pair because Fk, as defined by 
equation Q above, is a polynomial in the entries of 
the matrix A in which all terms are of order or 1 in 
any given matrix entry Ay (that is, F k is a "multi- 
affine" function of the entries of A) . 

Finally, we are interested also in the shift 5X in the 
eigenvalues A under the perturbation Sc. Noting that 
the eigenvalues A must be roots of the characteristic 
polynomial, 

J2{Fn-k + SF n - k )(\ + 8X) k = (16) 



may be used to derive 



SA = 



-£<SF„_ fe A fc 
fc 

EfcF^A*- 1 

k 



(17) 



A similar equation could be derived from equation 



( 16 1 for the mixed term AS A, but in the interest of 
brevity, we will not consider it explicitly. Instead, 
we will note only that the functional dependence of 
the eigenvalues A on the characteristic polynomial im- 
plies that the variation ASX must be a function of the 
variations SF k , AF k , and A5F k (as well as the unper- 
turbed eigenvalue A). 

5 Sparse Matrices 

Consider a distribution of systems for which the 
structure G is a random digraph in which each arc 
is included with independent probability p a rc, and 
for which all included arcs (j i) G Eq have 
linearization weight Ay drawn from a given proba- 
bility density function p w t{x) independently of the 
weights of all other entries (arcs). Then with prob- 
ability (1 — p^rc) tiie weight wk of any particular k- 
permutation term (henceforward fc-term) K in equa- 
tion Q is 0, while with probability p\ TC , it is distri- 
bution according to the density function 



Pfc-wt (w K | K c E G ) 



(18) 



Pwt 



7c- 1 



i=l 



WK 



'fc-1 

n 



Pwt{Xi) 



d k ~h 
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since all k arcs making up K are independent and 
identically distributed. We further take p w t symmet- 
ric about (so that the mean of each matrix entry is 
0). Thus, we consider the mean-square expectation 
value for wk for any fc-term K: 



.rare 



v k 

rare 



w t+,k) 



(19) 



w 2 K p k -wt (wk I K c E G ) dw K 



,J ^ + fe) ) is defined with respect to the dis- 
tribution assuming all arcs in the fc-term considered 
are present in the matrix A. Note particularly that 



( + fc j i does not depend on p arc . Finally, we also 

note that any two fc-terms K\ and K.% are distinct if 
and only if each of K\ and K-i contains at least one arc 
not contained in the other: but in this case they must 
have zero covariance ((wk 1 wk 2 )) = 0: since the (zero- 
mean) weight of each arc is independent of that of all 
others. Noting also that (wk) = for all K £ Qk - 
and hence (Fk) =0 — the mean-square value of the 
fc th coefficient of the characteristic polynomial is 



<^ 2 > = E «) - 



Kee k 



^ ^ Pare 



J (+,k) 



k\p k a 



(20) 



Now consider adding an additional SA qr to the 
matrix entry A qr . This will result in a change SFk 
proportional to the sum of the weights of all k- 
terms K containing the arc (r — > q). As there are 
Nk-qr = { 7 f c Z'2+s rq ) ~ 1)' possible fc-terms passing 
through the arc (r — > q), each having probability 
p^ 1 of being present (given the presence of the arc 
(r — > q)), the mean-square expectation value of SFk 
is given by 



N, 



k;qrP. 



fc-1 



(+,fc-l) 



\SA n 



(21) 



(Note the covariance of 5Fk with Fk vanishes owing 
to the symmetry of p w t ■ ) Now consider adding addi- 
tional weight 8A qr and AA UV to the two arcs (r — > q) 
and (v — > u). Similar considerations then lead to 



ASF k ] 



(v—}u),(r— >q) 



(22) 



Nl. v k ~ 2 

JV k;qruvP aTC 



W 



( + ,k-2) 



[AA 



where the final term corresponds to the Nk- qr uv — 

(i-^)(i-^)(r 4 4 :t^t^fc^:)^- 2 ) !possible 

fc-terms passing through both (r — > q) and (v — > u). 



From equations ( 19 1-( 22 ) , we can thus deduce that 



dA L 



^ oc J {k - 1] while oc J (k - 2! 

OC Pare , wmie — — — (X p ar c 

OAq r OA uv 



(23) 

from which we can conclude that for denser matrices 
(higher values of p a rc), the ratio of second- to first- 
derivatives of the root-mean-square expectation value 
of the coefficients of the characteristic polynomial Fk 
with respect to arc weights is generally lower than for 
sparser matrices. 

Noting that under the variation A H> A + SA + A A, 



[ASF k ] 



\SF k ] 



E 

i,j,q,r,s,u,v,w 



(24) 



a 2 F k 



d 2 F k 



dAijdA qr dA su dA vm 

AAijAA su ((SA qr SA VW )) S 



E dAij dAg r ((SAijSAqr)) 

we see that decreasing the relative sizes of the sec- 
ond derivatives of F k compared to the first deriva- 
tives might be expected to result in decreased values 
of the negative term in equation (11 ) (with U = Fk) 
for the correlation Corr,5 (SFk + ASF)., SFk). 

At this point it should be noted that the condi- 



tions under which equations ( 18 )- ( 24 1 were derived 



are quite restrictive and unrealistic if taken to de- 
scribe the linearization matrices A associated with 
systems f modeling biological systems. The simple 
random digraphs taken as the network structures G 
(especially the treatment of diagonal one-loop arcs 
(i — > i) symmetrically with all other arcs (r — > q)) 
and the assumption of independence of arc weights 
both neglect important features of real dynamic sys- 
tems. In particular, under these assumptions it is al- 
most certain that the matrix A will have some eigen- 
values with positive real part, and hence will not de- 
scribe the linearization of a system f about a stable 
steady state. 

We do not attempt to address these points ana- 
lytically, but instead shift our attention to numerical 
investigation of specific systems with more complex 
(and realistic) features. This also allows us to con- 
sider another dynamic system property — the dis- 
tance Re(Ai s ) from the set of eigenvalues of A to the 
imaginary axis in the complex plane (here designated 
the eigenvalue stability) — which is of more immedi- 
ate use in determining system stability. 
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6 Arc-Derivative Ratios 
Stabilized Matrices 



for where 



We first consider the impact of requiring system sta- 
bility — i.e., all eigenvalues of the linearization ma- 
trix A must be in the left half of the complex plane 
— on the ratio of second derivatives of the determi- 
nant or the eigenvalue stability to the first derivatives 
with respect to the arc weights Aij. Here we enforce 
this requirement on each matrix A by subtracting 
(Re(Ai s ) + l)/ from A, where Ai s is (one of) the largest 
real-part eigenvalue(s) of A, so that the new version 
of A will always have eigenvalue stability -1. (From 
the point of view of the network G corresponding to 
the matrix A, this operation consists of adding all 
one-loop arcs to G with weight — (Re(Ai s ) + 1).) 

In order to characterize the magnitudes of second- 
relative to first-derivatives of both the determinants 
det(.A) = (— l)"^ 1 ^ and the eigenvalue stabilities 
Re(Ai s ) with respect to arc weights A^ for a given 
matrix A with randomly constructed network struc- 
ture G, we numerically evaluated the 25 first deriva- 
tives and 625 second derivatives of both metrics with 
respect to a set E pt b containing 25 arcs randomly cho- 
sen from the arc set Eq- The root- mean-square (rms) 
average of the 625 second derivative terms was then 
divided by the rms average of the 25 first derivative 
terms to yield a single value 




E 

(J -h), (r->g)£-E p tt> 



d 2 U 
dA zi dA a 



(25) 



(j-H)e-Eptb 



dU 
OA,, 



of this ratio for each matrix A for each of the two 
metrics (determinant and eigenvalue stability). 

As mentioned in section [4] above, when comparing 
the response of systems with widely varying struc- 
tures to parametric perturbation, a measure of the 
perturbation size is required. For our purposes, it is 
natural to use the induced change to the eigenvalue 
stability as an indicator of the relative magnitude of 
a given perturbation, especially given that we have 
already normalized all of our stabilized systems to 
have the same base eigenvalue stability value of -1. 
In the discussion to follow, we thus also present re- 
sults for (rms averages of) derivatives with respect to 
the scaled arc weights 



A^ — 



ARe(A b 



Ai 



(26) 



ARe(Ai s ) 



\E 



ptb | 



>i)e-E ptb 



dRe(A ls ) 
dA in 



. (27) 

is approximately equal to the mean-square variation 
of Re(Ai s ) if the weights of the arcs in E pt b are 
subjected to independent Gaussian perturbations of 

2 

equal (small) variance ^ . The resulting scaled 
derivatives are then 



,>r -([ARc(A ls )f 



OA 



dU 

din 



(28) 



and 



d 2 U 1 / rAO a \-i d 2 U 

ajLel C = ?\ [ARe(Als)] ) Eptb dA-M, 



(29) 



where U is either the determinant ±_F n or the eigen- 
value stability Re(Ai s ). This scaling is specifically 
indicated where it has been applied. Clearly, for 
U = Re(Ai s ), one result of such scaling will be to 
bring the rms average of the first derivatives to unity. 

Figure[3]then shows the first, second (median), and 
third quartiles of these ratios calculated from popula- 
tions of stabilized random digraphic-structured ma- 
trices with varying arc density (100 matrices, each 
with a different structure G and n = 100, were gen- 
erated for each arc density grouping). Both the un- 
sealed and scaled data exhibit the expected trend of 
decreasing rms magnitudes of second derivatives com- 
pared to first derivatives for metric U chosen to be 
either determinant or eigenvalue stability. 

While these results suggest that that the qualita- 
tive effects of network density suggested by equa- 



tions (23) still hold for stabilized random digraphi- 



cal systems with varying arc densities (and perturba- 
tion sensitivities), we have not yet considered more 
complex structural variation. We also have not yet 
considered correlation under variation of the offset 
parameters originally introduced in section |4j as op- 
posed to the derivatives with respect to individual 
arc weights considered in this section. The behavior 
of these correlations (for random digraphical systems 
of varying arc density, for scale-free systems charac- 
terized by large node degree heterogeneity, and for 
highly clustered small-world systems) is investigated 
in the next section. 
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Figure 3: Median values of the ratio described by equation \25\ from randomly generated matrices (100 at each point) with 
random digraph structure of varying average node degree. Ratios are presented unsealed or scaled according to equations 
(|28|l-|29|l, as indicated. Y-axis labels for unsealed ratios on left, scaled to right. Error bars represent first and third quartiles. 



7 Correlation Susceptibility in 
Perturbed Dynamic Systems 

In order to study the effects of network structure on 
signal integration in non-linear dynamic systems de- 
scribed by equation it is necessary to consider 
not only the network structure G and the lineariza- 
tion matrix A, but also the quadratic terms B. For 
the purposes of this paper, random systems were con- 
structed by first generating a random network G (ac- 
cording to one of the random digraph generation al- 
gorithms described in appendix |A|, followed by pro- 
gressively building up A and B as follows. 

For each arc (j — > i) G Eq, we generate a random 



"rate constant" 



from a Gaussian distribution of 



vanishing mean and unit variance independent of all 
other system parameters and reset 



Aij Aij + 

Biji 1 y Biji -\- T^ij 
Baj 1 y Buj -\- T^ij 



(30) 



The contribution of arc (j — > i) to the system de- 
scribed by equation (|30| can be taken to represent a 



reaction 



sp 4 + sp -> [1 + sign(r li )] sp, + sp 



(31) 



with mass-action kinetics (offset by a reaction sp^ — » 
[1 — sign(rjj)] sp 4 with equal magnitude rate con- 
stant), where Xi measures the standardized deviation 
of the population of species i from it steady state 
value. 

After all arcs in Eq have been accounted for, the 
system S is stabilized by adding a multiple of the 
identity matrix to A so as to bring Re(Ai s ) to -1, as 
described in section [6] above (thereby including all 
one-loops (i — > i) in the final structure G). 

We then investigate the effects of perturbations Sc 
and Ac, as defined by equation Q, on the determi- 
nant ±F n and eigenvalue stability Re(Ai s ). Taking 
Ac to be fixed and 8c to be a random vector with co- 
variance matrix a 2 1 (and vanishing mean), and defin- 
ing for any function U of the parameters (^4, B, c) of 
the system about its steady state, 



dU 



and 



w = d 2 U 
13 deidej 



d 2 U 
dej dci 



(32) 



(33) 
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the variances and covariances of the changes SU and 
ASU to the arbitrary function U resulting from the 
random perturbation 6c in the presence of the fixed 
perturbation Ac may be derived from 

f)TT flTJ 

((6USU)) s =J2^g^((^S Cj )) s (34) 
id 1 3 
= a 2 u ■ u 

(^A^)) 5 = .^ Ac, 

hj 

= a 2 u T WAc 

({A6U ASU)) S = * 2 Y: PL AcjAck 



i,3,k 



dcjdci dckdci 



= a 2 (WAc) ■ [WAc) 

As discussed in section [I] above, we use the corre- 
lation equation (11) as a proxy metric to investigate 



the propensity of a system towards signal integration. 
For the reasons discussed in section [6j we consider 
eigenvalue stability-normalized perturbations satisfy- 



ing 



[5Re(A ls )] 2 \ = 



(35) 



where is a constant across all systems considered 
Taking U = Re(A 
obtain 



;<5Re(A ls )] : 



and employing equation (34 1, 
<9Re (Ai s 



- 2 E 



dc 



thus implying 
" 2 =fe 



dRe(A ls ) 



(36) 



(37) 



where u cs is the gradient vector of the eigenvalue sta- 
bility U cs — Re(Ai s ) with respect to the dynamic off- 
set vector c. Similar normalization of the interacting 
perturbation Ac then requires: 



<(A Ci Ac,)) A = 



<l> 2 



(38) 



where $ is, again, a system-independent constant, 
and we now consider Ac to be drawn from a proba- 
bility distribution with vanishing mean (Ac) = and 
covariance matrix proportional to the identity (simi- 



lar to 5c above) . Equations ( 34 ) - ( |38| together with 
equation ( |11[ ) imply 

Corra (SU, SU + ASU) = 1 - <^ 2 Tu, A + H.O.T. (39) 



where the A-correlation susceptibility Tu-A is given 

by 



U;A 



1 

2$2 



(WAc) • (WAc) 
u • u 



u - WAc 



u • u 



(40) 

The quantity <1> 2 Y(/ : a represents the degree to which 
Cons(SU, SU + ASU) is reduced in the context of the 
interacting perturbation Ac normalized to produce 
a first-order change ARe(Ai s ) of mean-square magni- 
tude 



[ARe (Ais)f 



= $ 2 



A 



(41) 



The A-correlation susceptibility Y[/ ; a depends not 
only on the particular system under consideration but 
also on the choice of perturbation Ac. As we are in- 
terested ultimately in characterizing the propensity 
of systems themselves toward signal integration or 
independence, we now further average over the inter- 
acting signal Ac to define the (average) correlation 
susceptibility, 



T„ = (T 



(42) 



T[7 (equation (42 )) is henceforward referred to as sim- 



ply the correlation susceptibility. 

Figure [4] provides boxplots of the correlation sus- 
ceptibilities Tjj for random digraph-, scale-free-, 
small- world-, and lattice-structured systems (see ap- 
pendix[A]for details regarding directed network struc- 
tures) generated as described above. ANOVA mod- 
els of the logged correlation susceptibilities against 
the six network groupings shown in figure [4] indi- 
cate significant effects (determinant -£5,594 = 9.68 
with p < 10~ 8 and R 2 — 0.075, eigenvalue stability 
^5,594 = 40.41 with p < 10~ 16 and R 2 = 0.25). 

As expected, figure [4] shows that the random 
digraph systems (rdg, rdglO, and rdg20) have 
progressively lower determinant- and eigenvalue 
stability-correlation susceptibilities as the arc density 
increases. We can now also compare the more com- 
plex scale-free (sf) and small-world (sw and latt) 
structured systems, however: both of these forms 
of structure result in lowered correlation susceptibil- 
ity (determinant and eigenvalue stability) when com- 
pared with random digraphical systems of the same 
arc density. In section [HJ we pursue structural met- 
rics on the basis of the arc density arguments made 
in section [5] designed to predict correlation suscepti- 
bility. 
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Figure 4: Boxplots (1.5 IQR whiskers) of determinant and eigenvalue stability correlation susceptibilities described by equation 
( |42| for varying types of network structure. Random digraphical systems of average node degree 6, 10, and 20 indicated by rdg, 
rdglO, and rdg20, respectively, while scale-free, small-world, and ring-lattice digraphs (all of average node degree 6) indicated 
by sf, sw, and latt. n = 100 nodes for all systems. 



8 Relevant Structural Metrics 

In sections [5] - [7| it has been argued that increasing 
network arc density leads to decreased values of our 
proxy for signal integration, correlation susceptibil- 
ity (both of the determinant and of the eigenvalue 
stability). One might expect, however, that not all 
arcs have an equal degree of influence on correlation 
susceptibility; both the weight Aji of an arc [i — > j) 
and the location of its termini i and j in the network 
suggest themselves as important factors to consider. 
Thus we seek a scalar metric derived from the lin- 
earization matrix A in the form of a sum of normal- 
ized arc weights, with the normalization factor hope- 
fully capturing some of the influence of arc locality. 

It is not immediately obvious how to go about con- 
structing such a normalization factor from the lin- 
earization matrix A. For the sake of clarity and par- 
simony, we begin by postulating that, since system 
stability is largely a function of the least-stable eigen- 
value^) Ai s of A, the eigenvector (s) associated with 
Ai s may offer a useful first indication to the relative 
influence of the various system nodes in questions of 
stability. 

Let v be (one of) the (generally complex) least- 
stable eigenvector(s) of the linearization matrix A, 



normalized so that v • v = 1. We here define the 
normalized arc density p A by 

PA = J2 l%lkl 2 N 2 (43) 

(i-yj)£E G 

For systems constructed as described in section [7J 
there will almost always be either exactly one (if X\ s 
is real) or two (if Ai s is complex) least stable eigen- 
vectors. When Ais is complex, the two least stable 
eigenvectors will be complex conjugates of each other; 
this implies that the value of pa as defined by equa- 
tion ( |43[ ) will be independent of which eigenvector is 
chosen. 

Figure [5] shows the association of the cor- 
relation susceptibility and the normalized 
arc density {R 2 (log(T dot ), \og(p A )) = 0.10, 
R 2 (log(T cs ),log(pA)) = 0.35). Considering only 
R 2 , the normalized arc density thus appears to 
explain more of the influence of network structure 
on correlation susceptibility than the ANOVA 
(section [7| on the network groupings themselves 
(this is not a completely fair comparison, however, 
since the normalized arc density includes more 
detailed information arising from the arc weights 
Aij). Additionally, the normalized arc density is 
capable of differentiating network structures with 
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Figure 5: Determinant and eigenvalue stability correlation stabilities (equation ( |42| ) versus normalized arc density (equation 
( |43| ) for varying types of network structure. Random digraphical systems of average node degree 6, 10, and 20 indicated by 
rdg, rdglO, and rdg20, respectively, while scale-free, small-world, and ring-lattice digraphs (all of average node degree 6) 
indicated by sf, sw, and latt. n = 100 nodes for all systems. 



the same "raw" arc density: consulting table [T] note 
that scale-free and, especially, small-world networks 
generally have larger normalized arc densities than 
random digraphs of the same average degree. 

Equation ( |43| considers only the magnitudes of the 
components of the least-stable eigenvector(s) v. Con- 
sidering equation ([7]), however, we see that the char- 
acteristic polynomial (and hence ultimately the eigen- 
values) of the linearization matrix A depend on sums 
of products of the weights of arcs making up permu- 
tation fc-terms. If, for example, the weights A\ 2 and 
A21 change as the result of a perturbation by <L4 12 
and 8A21, the resulting change to the weight of the 



fc-term K = (1 «-> 2) is 

6w K = A 12 SA 21 + SA 12 A 21 (44) 

This quantity depends on the relative phases of the 
arc weight perturbations, which will in turn depend 
on the relative phases of the perturbations to the 
components of the system steady state vector. Thus, 
we might suspect that the phases, as well as the mag- 
nitudes, of the components of the least-stable eigen- 
vector^), contain information relevant to the corre- 
lation susceptibility of a system. 

Thus, letting 5 = [y4~ 1 ( J 4 -1 ) T ], define the eigen- 





Tdet 


T 

± OS 


Pa 


9 A 


rdg 


3.6 ± 3.1 


17 ± 19 


0.22 ± 0.05 


0.015 ± 0.010 


rdglO 


3.5 ± 3.0 


11 ± 13 


0.23 ± 0.05 


0.012 ± 0.009 


rdg20 


2.2 ± 2.4 


5.7 ± 7.6 


0.32 ± 0.05 


0.012 ± 0.008 


sf 


1.5 ± 1.7 


2.6 ± 3.6 


0.29 ± 0.08 


0.018 ± 0.016 


sw 


3.0 ± 3.0 


7.0 ± 10 


0.40 ± 0.21 


0.016 ± 0.013 


latt 


1.3 ± 1.1 


0.13 ± 0.20 


1.2 ± 0.39 


0.030 ± 0.026 



Table 1: Median ± MAD values for determinant and eigenvalue stability correlation susceptibilities (equation j42[ l), normalized 
arc densities (equation (|43|), and eigenvector contractions (equation (|45|) for the systems used to generate figuresHl-pl 
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Figure 6: Determinant and eigenvalue stability correlation susceptibilities (equation i ]4'2[ ) ) versus eigenvector contraction (equa- 
tion l |45| ) for varying types of network structure. Random digraphical systems of average node degree 6, 10, and 20 indicated 
by rdg, rdglO, and rdg20, respectively, while scale-free, small-world, and ring-lattice digraphs (all of average node degree 6) 
indicated by sf, sw, and latt. n = 100 nodes for all systems. 



vector contraction by 



*A = 



Tr(H) 



(45) 



(Equation ( 45 1 is again generally independent of the 



choice of least-stable eigenvector when Ai s is complex, 
since the absolute value is invariant with respect to 
complex conjugation of its argument.) Note that ^ a 
is a weighted sum of the (phased) components of the 
least stable eigenvector v, with the weightings pro- 
vided by the diagonal elements of the matrix S. This 
choice of weighting was motivated by the observation 
that 



(46) 



(derived from Cov(i5c) = a 2 1 and equation Q). That 
is, we weight more heavily those components of the 
least stable eigenvector corresponding to nodes whose 
components in the steady state vector y undergo 
larger variance under the perturbations 5c. 

Figure [6] plots the correlation susceptibil- 
ity against the eigenvector contraction metric 
(i? 2 (log(T dot ),^) = 0.19, i? 2 (log(T es ),* A ) = 



0.27). Once again, the R 2 values of this metric are 
larger than those of the ANOVA from section [7] on 
the network groupings — especially, in this case, 
with regard to the determinant correlation suscep- 
tibility. Compared to the normalized arc density, 
however, the eigenvector contraction has a relatively 
large within-network-group MAD-to-median ratio, 
indicating that this factor perhaps helps to explain 
more of the within-group variance in correlation 
susceptibility. 

The correlation of the eigenvector contraction and 
the two correlation susceptibilities is negative, simi- 
lar to the relationship between normalized arc den- 
sity and correlation susceptibilities. Following the 
arguments made in section [5] one might consider 
the relative impact of in-phase- versus out-of-phase- 
perturbations to the arc weights (originating from 
perturbations Sy to the system steady state) on the 
first- as compared to the second-derivatives of the 
fc-terms composing the characteristic polynomial Fk- 
For the sake of brevity, this paper does not continue 
on this path, noting only the empirical correlation 
between contraction and susceptibility. 

Table [2] provides the results of regression fits 

log(f ) = ff oe(pA j logG^ ) + P%J>A (47) 
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Figure 7: Determinant and eigenvalue stability correlation susceptibilities (equation | |40[ l) versus regression value (equation 
XXX) for varying types of network structure. Random digraphical systems of average node degree 6, 10, and 20 indicated 
by rdg, rdglO, and rdg20, respectively, while scale-free, small-world, and ring-lattice digraphs (all of average node degree 6) 
indicated by sf, sw, and latt. n = 100 nodes for all systems. 





Estimate 


Std Error 


Pr(> \t\) 


adet 

P log(p A ) 
odct 
P* A 


-0.171 
-0.363 


0.040 
0.040 


1.87E-5 
6.66E-19 


acs 
fit* 


-0.459 


0.033 


3.53E-38 


-0.339 


0.033 


7.67E-23 



Table 2: Standardized regression coefficients for models de- 
scribed by equation l |47| . The resulting fits have R^ Et = 0.21 
and Rl =0.45. 



for U cither determinant or eigenvalue stability; here 



we use the tilde notation X 



x- 



for a random 



variable X with estimated mean fix and estimated 
standard deviation &x ■ Figure [7] plots (on the orig- 
inal scale) the true correlation susceptibility values 
against their regression values. 

9 Conclusions 

This paper considers the behavior of the correlation 
of two quantities, the determinant det(A) and the 
eigenvalue stability Re(Ai s ), associated with the sys- 
tem linearization A under the influence of perturba- 
tions of the system dynamics of the form equation 



These correlations are presented as indicators 
of the degree to which distinct signal inputs inter- 
act cooperatively versus acting independently. The 
necessity of normalization of the perturbation mag- 
nitudes leads to the introduction of the correlation 



susceptibility (equations (39 1 - (42 1). 



Sections [5] - [6] suggest that increasing network arc 
density may decrease the determinant and eigenvalue 
correlation susceptibilities. Section [7] confirms this 
suggestion via simulation for a class of quadratic dy- 
namic systems, while also indicating similar trends 
with regard to network clustering and degree hetero- 
geneity. Section [8] constructs a structural metric, the 
normalized arc density (equation (43)), which pro- 



vides a potential unified explanation for the impacts 
of these distinct aspects of network structure. 

The results with regard to (normalized) arc den- 
sity are somewhat reminiscent of the dimensionality 
dependence of mean field theory in statistical physics 
7\. An Ising model spin is less sensitive to fluctua- 
tions of a nearby spin in a higher dimensional/higher 
connectivity lattice in a manner similar to the de- 
creasing sensitivity of the response to one perturba- 
tion with respect to the influence of another in the 
dynamic systems studied here as normalized arc den- 
sity increases. 
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Section [8] also introduces the eigenvalue contrac- 
tion metric, which appears to offer a predictor of the 
correlation susceptibilities complementary to the nor- 
malized arc density. Both the network density and 
the eigenvalue contraction focus on the nodes most 
involved in the least stable eigenvector(s) of the lin- 
earization A, which may be loosely thought of as the 
points at which the network is most susceptible to 
destabilization. 

The results with regard to the normalized arc den- 
sity and eigenvector contraction indicate that net- 
works with lower arc density in the neighborhood of 
these focal nodes, and for which the fluctuations of 
the dynamic variables associated with these nodes are 
anticorrelated/out-of-phase with regard to the least 
stable dynamic modes, are more likely to respond to 
multiple signal inputs in a cooperative manner. It 
is intriguing to contemplate extension of this idea to 
develop more easily calculable (i.e., not linearization 
A-eigendccomposition-dependent) system structural 
metrics for potential use as predictors of the degree 
of interaction between distinct perturbations (such 
as, for example, an infection and a drug treatment). 



Appendices 



A Generation of Directed Net- 
work Structures 

Random digraph structures were generated by inde- 
pendently adding each arc (i — ¥ j) to the digraph with 
a fixed probability p mc = , where the parameter 
d controls the average node in- and out-degree and n 
is the number of nodes in the network. 

Directed scale-free networks with n nodes and 
2nd arcs were generated in a manner similar to 
the Barabasi preferential attachment mechanism [5]. 
First, (d+1) fully connected nodes were added. Then 
n— (d+1) nodes were added sequentially one-by-one, 
with d arcs added directed from the newly added node 
to old nodes and d arcs added from old nodes to the 
new node at each step. For each of these 2d new arcs, 
the identity of the adjacent old node (whether it be 
tail or head of the new arc) was chosen at random 
with non-uniform probability proportional to the sum 
of the in- and out-degrees of the old node (subject to 
the constraint that no two arcs may connect the same 
pair of nodes with the same directionality) . 

Small-world digraphs with n nodes and 2nd arcs 



(where d is even) were generated starting with a 
ring lattice in which each node is connected bi- 
directionally to the | nearest nodes in each direction, 
then rewiring each directed arc with probability p IV! . 
Similar to the mechanism described by Watts, et. al. 
..9], each rewired arc has one terminus (head or tail 
chosen randomly with uniform probability) changed 
to a node chosen randomly from all nodes in the sys- 
tem with uniform probability, again constrained to 
prevent any two arcs from sharing both the same tail 
and the same head. 
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